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Abstract 

Lattice Boltzmann simulations are used to explore the behavior of liquid 
crystals subject to Poiseuille flow. In the nematic regime at low shear rates 
we find two possible steady state configurations of the director field. The 
selected state depends on both the shear rate and the history of the sample. 
For both director configurations there is clear evidence of shear-thinning, a 
decrease in the viscosity with increasing shear rate. Moreover, at very high 
shear rates or when the order parameter is large, the system transforms to a 
log-rolling state with boundary layers that may exhibit oscillatory behavior. 

I. INTRODUCTION 

As suggested by their name, liquid crystals are fundamentally a liquid and hence hydro- 
dynamics can be very important to their behavior. Moreover, because they comprise rod-like 
molecules the properties of their flow can be much richer than that seen in simple fluids. 
This is because the translational motion of the fluid is coupled to the inner orientational 
motion of the molecules. As a result the flow disturbs the alignment. Conversely changing 
the alignment will almost invariably induce a flow. 

This coupling can have important practical consequences in applications of liquid crys- 
tals. Back-flow can be an important limiting factor in determining the rate at which a 
liquid crystal display device can be switched [0 . The ordering induced by a shear flow could 



potentially be beneficial if understood and controlled or can be the origin of mechanical 
instabilities under processing conditions @||. 

In addition to technologically relevant issues, there are also a number of fundamental 
scientific questions which can be addressed in the context of liquid crystals. It is known that 
shear flow can induce a phase transition from the isotropic (disordered) to nematic (liquid 
crystalline ordered) phase The resulting non-equilibrium phase diagram includes a line 

of first-order transitions ending in a critical point. Techniques for examining and classifying 
such non-equilibrium phase transitions are still in need of further testing and development. 
In addition, a large number of structured fluids, like a liquid crystal, have order parameters 
that have tensor characteristics. Methods for studying such systems are still somewhat 
lacking and need to be explored. 

Although the effects of hydrodynamics are well known, it is usually very difficult to 
fully incorporate these effects into a calculation. This has been due to the complexity of 
the nine coupled non-linear partial differential equations which describe the system (the 
order parameter is a symmetric traceless tensor with five components, the Navier-Stokes 
equation has three components, and the continuity equation is scalar). While the final 
results may appear simple it can often be an imposing task to derive them analytically from 
the underlying equations. From a numerical viewpoint accurate simulations can be difficult 
as the solutions can contain topological defects requiring resolution of length scales much 
smaller than those associated with hydrodynamic flows. 

In an attempt to overcome these problems we have developed a modified lattice Boltz- 
mann algorithm to simulate liquid crystals. The lattice Boltzmann algorithm, which includes 
hydrodynamics, makes use of a physical analogue (from statistical mechanics) to map the 
partial differential equations onto equations for the evolution of two probability distribu- 
tions. Moments of these distributions are then related to the physical variables of interest. 
Among the advantages of this reformulation is that the resulting equations are much sim- 
pler to model and possess greater numerical stability than, say, a traditional finite-difference 
scheme applied to the original equations. 
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In this paper we apply the lattice Boltzmann method to the study of Poiseuille flow in 
liquid crystals. The paper is organized as follows. First, we describe the model of liquid 
crystal hydrodynamics proposed by Beris and Edwards . This contains both the Erickson- 
Leslie and the Doi models as limiting cases. Second, we outline the main features of the 
lattice Boltzmann algorithm used to simulate the model. Results are presented to show 
the effect of the flow on the director orientation and how this in turn leads to both shear 
thinning and hysteresis. We also examine the transition to the log-rolling state. 

II. THE EQUATIONS OF MOTION FOR LIQUID CRYSTAL HYDRODYNAMICS 

We follow the formulation of liquid crystal hydrodynamics described by Beris and Ed- 
wards |7j . The continuum equations of motion are written in terms of a tensor order param- 
eter Q. The advantage of this approach is that it includes both the isotropic (Q = 0) and 
the nematic (Q 7^ 0) phases and allows an order parameter of variable magnitude within the 
latter. Hence it is possible to explore the effect of flow on the phase transition between the 
two states. Moreover the hydrodynamics of topological defects (point defects in two dimen- 
sions) are naturally included in the equations. The Beris-Edwards equations reduce to the 
Ericksen-Leslie formulation || of nematodynamics for a uniaxial nematic in the absence of 
defects. 

We describe the equilibrium properties of the liquid crystal by the Landau-de Gennes 
free energy || 

T = \d\ l~Ql p - h -Q a pQ M Q, a + |(Q^) 2 + 1{d a Qp X f j , (II. 1) 

where Greek indices represent Cartesian directions and a sum over repeated indices is as- 
sumed. This free energy describes a first-order transition from the isotropic to the nematic 
phase. Note that, for simplicity, we are working within the one elastic constant approxima- 
tion. Three elastic constants are needed to fully characterize the nematic phase ||. 
The equation of motion for the nematic order parameter is 

(d t + u- V)Q - S(W, Q) = EH (II.2) 
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where T is a collective rotational diffusion constant. The first term on the left-hand side of 
equation ( |11.2| ) is the material derivative describing the usual time dependence of a quantity 
advected by a fluid with velocity u. This is generalized by a second term 

S(W, Q) = (£D + n)(Q + 1/3) + (Q + 1/3) (^D - ft) 

-2f (Q + I/3)Tr(QW) (II.3) 

where D = (W + W T )/2 and f2 = (W — W T )/2 are the symmetric part and the anti- 
symmetric part respectively of the velocity gradient tensor W a p = d$u a . S(W, Q) appears 
in the equation of motion because the rod-like shape of the molecules means that the order 
parameter distribution can be both rotated and stretched by flow gradients. £ is a constant 
which depends on the molecular details of a given liquid crystal. 

The term on the right-hand side of equation |II.2j describes the relaxation of the order 
parameter towards the minimum of the free energy. The molecular field H which provides 
the driving force is related to the derivative of the free energy by 

H -{%- mTr %] (IL4) 

= -aQ + b (Q 2 - (I/3)TrQ 2 ) - cQTrQ 2 + kV 2 Q. 
The fluid, of density p obeys the continuity 

d t p + d a pu a = (II. 5) 

and the Navier-Stokes equations 

pd t u a + pupdpu a = dpT a p + dpaap (H-6) 
+^-(df3((S a/3 - 39 p P )5 7 m 7 + d a up + dpUa). 

where Tf is related to the viscosity and P is the pressure, 

^o = pT-|(VQ) 2 . (II.7) 
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The form of this equation is not dissimilar to that for a simple fluid. However the details of 
the stress tensor reflect the additional complications of liquid crystal hydrodynamics. There 
is a symmetric contribution 

<y a /3 = —PoS a /3 — iH (X1 {Q 1 p + — 5 7( g) 

1 1 

(II.8) 

and an antisymmetric contribution 

T af3 = QayHyp — H al Q 1 p. (H-9) 

A detailed account of the theories of liquid crystal hydrodynamics can be found in the book 
by Beris and Edwards J7J. 

III. A LATTICE BOLTZMANN ALGORITHM FOR LIQUID CRYSTAL 

HYDRODYNAMICS 

We now define a lattice Boltzmann algorithm which solves the hydrodynamic equations of 
motion of a liquid crystal ( |11.2|) , ( [11. 5|) , and ( |H.6|) . Lattice Boltzmann algorithms are defined 
in terms of a set of continuous variables, usefully termed partial distribution functions, 
which move on a lattice in discrete space and time. For a simple fluid a single set of partial 
distribution functions which sum on each site to give the density is needed. For liquid crystal 
hydrodynamics this must be supplemented by a second set, which are tensor variables, and 
which are related to the tensor order parameter Q. 

We define two distribution functions, the scalars fi(x) and the symmetric traceless 
tensors Gj(x) on each lattice site x. Each fi, Gj is associated with a lattice vector 
e*j. We choose a nine-velocity model on a square lattice with velocity vectors e*j = 
(±1, 0), (0, ±1), (±y/2, ±y/2), (0, 0). Physical variables are defined as moments of the distri- 
bution function 

P = Y,fi' pUa = J2fi e ic" Q = G * (HI.10) 

i i i 
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The distribution functions evolve in a time step At according to 

fi(x + eiAt,t + At)- fi(x,t) = (111.11) 
At 

— [C fi {S, t, {fi}) + C fi {S + e t At, t + At, {/*})] , 

Gi(x + eiAt,t + At) -Gi(x,t) = (III. 12) 

At 

— [C Gi {S, t, {d}) + C Gi (x + e { At, t + At, {G*})] 

The left-hand side of these equations represents free streaming with velocity e*j, while the 
right-hand side is a collision step which allows the distribution to relax towards equilibrium. 
/* and G* are first order approximations to MS + e^Ai, t + At) and Gj(x + e*jAt, t + At) 
respectively. They are obtained from equations ( |111.11| ) and ( [111.12j ) but with /* and G* 
set to fi and G». Discretizing in this way, which is similar to a predictor-corrector scheme, 
has the advantages that lattice viscosity terms are eliminated to second order and that the 
stability of the scheme is improved. 

The collision operators are taken to have the form of a single relaxation time Boltzmann 
equation, together with a forcing term, 

C fl (x,t,{f l })= (111.13) 

--(Ms, t) - n\s, t, {fi})) + Vi {s, t, {f t }), 

C Gl (x,t,{G t })= (111.14) 
-— (Gi(S, t) - Gf(x, t, {G,})) + Mi{S, t, {Gi}). 

The form of the equations of motion and of thermodynamic equilibrium follow from the 
choice of the moments of the equilibrium distributions fl q and G^ g and the driving terms 
Pi and Mj. fl q is constrained by 

= Pi H fi 9e i<* = pU a , fi 9e iaCiP = -Va/3 + pU a U{5 (III. 15) 

i i i 

where the zeroth and first moments are chosen to impose conservation of mass and momen- 
tum. The second moment of f eq controls the symmetric part of the stress tensor, whereas 
the moments of Pi 
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E^ = > ^PiZia = dpTap, ^2pie ia eip = (III. 16) 

i i i 

impose the antisymmetric part of the stress tensor. For the equilibrium of the order param- 
eter distribution we choose 

E G i Q = Q> E Gi q e ia = Qu a , G- 9 e ia e i/3 = Qu a up. (III. 17) 

i i i 

This ensures that the order parameter is convected with the flow. Finally the evolution of 
the order parameter towards the minimum of the free energy is most conveniently modeled 
by choosing 

J2 M 4 = TH(Q) + S(W, Q), J2 M * e - = (E M *K- ( IIL18 ) 



Conditions ( Jill. 15[ ) ( [L11.18| ) can be satisfied as is usual in lattice Boltzmann schemes by 



writing the equilibrium distribution functions and forcing terms as polynomial expansions 
in the velocity 0. 

Taking the continuum limit of equations (|III.11|) and (|III.12|) and performing a Chapman- 



Enskog expansion leads to the equations of motion of liquid crystal hydrodynamics ( |II.2| ), 
flini) , and ( TO . 



IV. POISEUILLE FLOW 

We now consider Poiseuille flow in a liquid crystal. Poiseuille flow is the flow between 
two non-slip boundaries at y = and y = h, driven by a pressure gradient or body force. In 
a simple fluid this geometry leads to a parabolic profile of velocities [ITU 



Ux = ~{pbf)v(v-h) (IV.19) 

where rj is the viscosity, (pbf) is an applied body force, and u x is the velocity along the 
channel. A simple example is gravity-driven flow. 
In an experiment the volumetric flow rate 

*u x dy = ^-(pb f )h 3 (IV.20) 
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can be measured as a function of the applied body force (or pressure gradient) to obtain 
the viscosity rj. If the material being studied is homogeneous, then the viscosity should be a 
unique function of the effective shear rate, Z/h 2 . This can be used as a working definition of 
a simple fluid. If the fluid is not homogeneous, then the relation between rj and the effective 
shear rate Zj h 2 will depend on geometrical factors such as the channel width h. 

In particular, a nematic liquid crystal has internal structure which can lead to inhomo- 
geneities. The internal structure is strongly affected by the presence of the flow field. In 
a simple shear flow the director in a nematic liquid crystal will try to align with the flow 
direction at an angle 6 given by 

£cos2£ = -^-, (IV.21) 
2 + q 

where q is the amplitude of the order parameter (equal to 3/2 times the magnitude of the 
largest eigenvalue of Q) QJ. However this condition is not compatible with fixed boundary 
conditions at the walls of the channel and therefore the direction of the director field will be 
determined by a balance of competing factors. In addition note that because £ < 1 there will 
be no solution to equation (|lV.21p as q approaches one. When no solution is available the 



director tumbles in the flow. An examination of the different regimes of tumbling under shear 
(but imposing a constant shear gradient across the sample rather than solving Navier-Stokes 



equations with shear boundary conditions) is given in Ref. 

In order to make it easier to compare our results to a wider array of experimental 
situations, we will plot all quantities in dimensionless numbers. First, note that 1/r has 
units of viscosity and k has units of force. One can then form a velocity scale as kY/H. 
Using this velocity scale one can construct the Erickson number Er = u x h/(nT) (normally, 
in the Erickson-Leslie theory, the Ericksen number is defined using a viscosity. However, it 
can be shown that V sets a scale for the contribution to the measured viscosity from the 
liquid crystal order |7[], which is what we are interested in here). The Ericksen number for 
the flows we examine is in the range 10 2 to 10 4 . For comparison, the Ericksen number for 



the flows in the experiments described in Reference |L3[ is in the range 10 to 10 . A useful 



time scale is l/(poT) where po is the ambient pressure (po — latm in all our simulations). 
All the simulations are performed in two spatial dimensions, while the order parameter is in 
three (i.e. the director can point out of the plane). 

The remainder of the paper is organized as follows. We first examine the order parameter 
configurations which are possible in different flow regimes. Then the corresponding effective 
viscosity is measured as a function of flow rate. We find shear thinning, and that the viscosity 
is not a unique function of the effective shear rate, in good agreement with experiments. 
Applying a scaling suggested by Ericksen |L2] we find that the data can be collapsed onto 



a single curve with two branches. The transition to the log-rolling state, where the director 
points out of the shear plane, is then briefly examined. 

A. Order parameter profile 

Distortions in the nematic state are induced by the laminar Poiseuille flow. Two possible 
configurations of the director field, for the case with strong normal anchoring at the bound- 
aries, are shown in Figure |J In ) the system was initialized in the nematic phase 
with the director uniform and perpendicular to the boundary before the flow was turned 
on. In case (b) the system was kept in the isotropic phase for the first 4000 time steps to 
establish the flow and then quenched into the nematic phase. Otherwise all parameters were 
identical. It is clear that the two systems have chosen different compatibility conditions at 
the center of the flow. Configuration (a) is described in the book by deGennes M . Case (b) 
does not appear to have been previously considered in the context of flow. It is, however, 
very similar to the equilibrium director configuration expected in a cylindrical tube with 
normal anchoring boundary conditions. 

To understand these configurations note that in both cases there are three distinct regions 
of the flow. First, away from the walls and the center, the director is oriented at the angle 
given by equation ( |IV.21|) . Approaching the wall, the director must change its configuration 
to match the boundary condition (normal anchoring). The distance e\ over which this 
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occurs is determined by a competition between the elastic energy and the shear stress at the 
boundary and can be estimated as 

ei ~ (K/t] Si ) 1/2 (IV.22) 

where K is an elastic constant, rj a viscosity and s\ is the shear rate at the boundary ||. In 
terms of the parameters of our model this becomes 

d ~ (kT/ Si ) 1/2 (IV.23) 

as the polymer viscosities can be shown to be proportional to 1/T (?]]. 

Similarly at the center of the flow the director must turn to extrapolate between the 
directions on the two sides. This will happen over a distance 

e 2 ~ {kT/s 2 ) 1/2 . (IV.24) 

where s 2 is the shear rate at a distance €2/2 from the center. 

For very low (or zero) flow rates the boundary and central regions will join and the 
configuration (a) with the director perpendicular to the flow at the center will be preferred 
(this should be clear for the zero flow case where this configuration has no elastic energy). 
At high flow rates or for wider channels the configuration with the director aligned with the 
flow at the center will be stable. To see why this is so, consider the elastic energy in the 
central region. Outside the central region, the director makes an angle ±8 with the flow 
direction. As can be seen from equation ( |IV.21| ), this angle is always less than 45° and it will 
therefore cost less energy to go through 0° in the center than 90°. Hence the flow aligned 
configuration (b) is preferred. Note, however, that (a) can exist as a metastable state, given 
suitable initial conditions. 

As the flow rate is increased the system can spontaneously change the director configu- 
rations at the center from perpendicular to parallel to the flow. This happens by nucleating 
a short lived defect at the center which "unzips" the director. Figure [2] shows the amplitude 
of the order parameter at the center of the flow as a function of time for a situation where 
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this happens. The sharp dip corresponds to the appearance of the defect. To correctly 
describe such a change of state clearly requires a formalism, such as that used here, where 
the amplitude of the order parameter can vary. 



B. Shear thinning 

The viscosity of a liquid crystal depends on its orientation with respect to the flow ||. 
As the orientation of the director at the boundaries and center of the channel depends on 
the flow rate the viscosity will be a function of the flow rate. This dependence is shown in 
Figure [| The apparent viscosity was calculated assuming equation( |1V.19| ) holds (indeed 
the velocity profile does remain close to parabolic). This is equivalent to the procedure used 
in experiments ]13| where equation ( |IV.20|) is assumed to hold and the apparent viscosity 



is measured SIS db function of flow rate. 

As expected the apparent viscosity also depends on the configuration of the director at 
the center of the flow and hence the viscosity curves in Figure 3 have two branches. As the 
flow and hence the strain rate is increased, the viscosity decreases: the fluid exhibits shear 
thinning. At high shear rates the curves approach each other as the boundary and central 
regions become smaller and the viscosity is dominated by the flow-aligned regions. 



Ericksen has shown ||12||, based on a dimensional argument, that the apparent viscosity 



at a shear rate s can be written 



Vapp(s) = T) app (0)f 



h 



Vapp(0)f 



kT\ 1/2 1 
Si J h 



(IV.25) 



where the dimensionless function / depends on the ratio between various Leslie coefficients 
and the particular laminar flow under study. Scaling the data in this way, as shown in 
Figure |], we find that the results from the two different channel width do indeed collapse 
onto a single curve, with two branches for the two different director configurations. 
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C. Log-rolling state 



So far we have assumed that the director stays in the shear plane. This is not always the 
case. If no solution to equation ( |1V.21| ) exists, the director tends to move out of the plane 



to a log-rolling state. In this state the director is stationary and perpendicular to the shear 
plane. However, the eigenvector corresponding to the second largest eigenvalue of Q (the 
sub-director), remains in the shear plane but oscillates (thus the "rolling" in log-rolling). As 
one approaches the boundary where the director is constrained to be in the shear plane (and 
perpendicular to the wall) the director must rotate back into the plane. In this regime, the 



director in the boundary layer may be stationary or may exhibit oscillatory states ||11|| . An 
example of one of these dynamic flows is shown in Figure [|. The oscillation of the director 
in the boundary layer is plotted explicitly in Figure || There remain interesting questions 
about the nature of the transitions between different types of director behavior in shear 
flows. 

V. DISCUSSION 

We have extended lattice Boltzmann algorithms of multiphase flow to treat the case of a 
non-conserved, tensor order parameter. Hence it is possible to simulate the Beris-Edwards 
equations of liquid crystal hydrodynamics. These hold for the isotropic, uniaxial nematic, 
and biaxial nematic phases. 

Using the approach to investigate Poiseuille flow in liquid crystals a rich phenomenology 
is apparent because of the coupling between the director and the flow fields. In weak shear, 
for strong normal anchoring at the boundaries, two states can be stabilized. For narrow tubes 
the director field at the center of the flow prefers to align perpendicular to the boundaries, 
for wider tubes it prefers to align along the flow direction. The effective viscosity decreases 
markedly with increases in the shear rate. This is in agreement with the shear thinning seen 



in experiments [13 



As the shear rate is increased the stable configuration changes to a log-rolling state where 
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the order parameter tumbles with the flow. Oscillations in the director orientation in the 
boundary layer are also observed. In experiments topological defects are nucleated at these 



shear rates, perhaps because of impurities p^JT5[] . In the simulations a similar phenomenon 
may occur if fluctuations are added and work is in progress in this direction. 

There are many avenues for further research opened up by the rich physics inherent in 
liquid crystal hydrodynamics and the generality of the Beris-Edwards equations. The pos- 
sibility of numerical investigations will prove very helpful as the complexity of the equations 
makes analytic progress difficult. Of particular interest are the possibilities of exploring non- 
equilibrium phase transitions such as shear banding: the coexistence of states with different 
strain rates. In particular it will be possible to investigate the pathways by which shear 
bands form. Work is also in progress to investigate the effects of flow on the dynamics of 
topological defects. 
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FIGURES 




FIG. 1. Left: Two different director configurations for steady-state Poiseuille flow with the 
director in the plane. The boundary conditions are such that the director is perpendicular (9 = 90) 
to the walls at the top and bottom of the figure and the flow is from left to right. For clarity, only 
the director on every third lattice point is displayed. Right: The corresponding fluid velocity v y . 
The solid and dotted lines are for the director configurations (a) and (b) respectively. Velocities 
are scaled by nT/h to make them dimensionless. 
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FIG. 2. The amplitude of the order parameter q at the center of the channel as a function 
of time t, measured in units of l/(poF), for a flow strong enough to cause the system to switch 
between the two states (a) and (b) in Figure 1. This occurs via a two-stage process. The distorted 
region at the center is very slowly compacted, then quickly 'unzipped' by a defect running along 
the channel. 
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FIG. 3. Apparent viscosity ijasa function of the average strain. (The solid curves correspond 
to channels twice as wide as the dotted curves.) Different branches correspond to the director 
perpendicular (upper branch) or parallel (lower branch) to the flow at the center of the channel. 
The strain rate is scaled by poT and the viscosity by 1/T to make them dimensionless. 
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FIG. 4. Apparent viscosity as a function of the scaled variable (hy/s)^ 1 . (k and the viscosity 
coefficients are held the same for the two cases, only the channel width is changed.) Different 
branches correspond to the director perpendicular (upper branch) or parallel (lower branch) to the 
flow at the center. The units are the same as for Figure 3. 
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FIG. 5. Steady and oscillating states in Poiseuille flow. The lines represent the director orienta- 
tion (eigenvector corresponding to largest eigenvalue of Q) projected down onto the xy-plane, and 
shading represents the amplitude of the order parameter (largest eigenvalue) . Flow is from top to 
bottom, and the walls are at the left and right. At the walls, the director is aligned perpendicular 
to the boundary, (a) Configuration (b) from Figure 1. (b) Snapshots of an oscillating configuration 
where the central region is in the "log-rolling state" (director perpendicular to the plane) and the 
boundary region consists of a transition from a configuration in the shear plane to a "tumbling" 
(director rotating in the shear plane) and "kayaking" region (director rotating in and out of the 
plane) interfacing to the central log-rolling state. 
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FIG. 6. Component of the director in the direction of fluid flow as a function of time at a point 
in the oscillating boundary layer (see Figure 5). 
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